/* 
	 Example of a hybrid API using OpenMP and MPI 
	 . Sums the elements of a 1D array of size n.
   . Uses openMP reduction and MPI_Reduce()
	 
	 using a gcc-based MPI
	 > mpic++ -fopenmp reduce.cc -o reduc
   using a icc-based MPI
	 > mpic++ -openmp reduce.cc -o reduc
	 and run with
	 > qsub -I -l nodes=2
	 > export OMP_NUM_THREADS=4; mpirun -np 2 ./reduce [array_size]
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>
#include <vector>

using namespace std;

#define N_DEFAULT 1000
#define REAL float
#define USE_GPU 0
#define MAIN_RETURN_NOERROR 0

int main( int argc, char *argv[]){

	int n; // problem size
	int namelen;
	char name[32];

	MPI_Init(&argc,&argv); 	// intialize MPI and MPI-related variables 
	MPI_Comm comm=MPI_COMM_WORLD;
	int rank; MPI_Comm_rank( comm, &rank);
	int size; MPI_Comm_size( comm, &size);
	MPI_Get_processor_name(name,  &namelen);
	int root = 0;
	bool iamroot = (rank==root);

	// read problem size and broadcast to all processors
	if(iamroot)	n= (argc==2) ? atoi(argv[1]) : N_DEFAULT;
	n = int(floor(double(n)/size));
	MPI_Bcast((void*)&n, 1, MPI_INT, root, comm);
	printf("Rank[%d], global size=%ld, local size=%d, host=%s\n", rank, long(n)*size, n, name);

	// initialize array values
	vector<REAL> A; 	A.resize(n);
#pragma omp parallel for
	for(int i=0;i<n;i++) A[i] =REAL(1e-2);  

	double tstart, tend;
	tstart = MPI_Wtime();
	// OMP reduction
	double sum;
#pragma omp parallel for reduction(+ : sum)	
		for(int i=0; i<n; i++) sum += double(A[i]);		

		// MPI Reduction
	double sumo;
	MPI_Reduce( (void *)&sum,  // input
							(void *)&sumo, // final result
							1,             // data length
							MPI_DOUBLE,    // data type
							MPI_SUM,       // reduction operator
							0,             // process id that will collect the answer
							comm           // communicator
							);
	tend = MPI_Wtime();

	int ntds;		
#pragma omp parallel
#pragma omp master
	ntds = omp_get_num_threads();

	//print  partial sum in process
	printf("Rank [%d], partial sum is: %.2g\n", rank, sum);

	if(iamroot) printf("\n\tUsing %d processes and %d threads/process\n", size,ntds);
	if(iamroot) printf("\tThe answer is %.2g. It should be:%.2g\n\n", sumo, double(n*size)*1e-2);
	
	MPI_Barrier(comm); 
	fflush(stdout);
	printf("Rank [%d], ellapsed time=%2.2g secs\n", rank, tend - tstart);
	fflush(stdout);
	MPI_Finalize();
	return MAIN_RETURN_NOERROR;
}




	
 

	
	
